clear
rng default
n_types=4;
all_years=1940:1:1967;
n_years=size(all_years,2);

%% Construct pictures reporting the bidimensional projections of the identified sets computed in the cluster. 
%When the hitting the grid bounds report -Inf, +Inf in the table below. 
%3 indicates -Inf, +Inf
%2 indicates lower bound -Inf
%1 indicates upper bound +Inf
%NaN indicates finite lower and upper bound

%U	2	3	4	2	3	4	2	3	4	2	3	4
U=[ 1	3	3	1	3	3	1	3	3	1	1	1; ...
 	1	3	3	1	3	3	1	3	3	1	1	1; ...
 	1	3	3	1	3	3	1	1	3	1	1	1; ...
 	1	3	3	1	3	3	1	1	3	1	1	1; ...
 	1	3	3	1	3	3	1	1	3	NaN	1	1; ...
 	1	3	3	1	3	3	1	1	3	NaN	1	1; ...
 	1	3	3	1	3	3	1	1	3	NaN	1	1; ...
 	1	3	3	1	3	3	1	1	3	NaN	1	1; ...
 	1	3	3	1	3	3	1	1	3	NaN	1	1; ...
 	1	3	3	1	3	3	1	1	3	NaN	1	1; ...
 	1	3	3	1	3	3	1	1	3	NaN	1	1; ...
 	1	3	3	1	3	3	1	1	3	NaN	1	1; ...
 	1	3	3	1	3	3	1	1	3	NaN	1	1; ...
 	1	3	3	1	3	3	1	1	3	NaN	1	1; ...
 	NaN	2	2	1	3	3	NaN	1	3	NaN	1	1; ...
 	NaN	2	2	1	3	3	NaN	1	3	NaN	1	1; ...
 	NaN	2	2	1	3	3	NaN	1	3	NaN	1	1; ...
 	NaN	2	2	1	3	3	NaN	1	3	NaN	1	1; ...
 	NaN	2	2	1	3	3	NaN	1	3	NaN	NaN	1; ...
 	NaN	2	2	1	3	3	NaN	1	3	NaN	NaN	1; ...
 	NaN	2	2	1	3	3	NaN	1	3	NaN	NaN	1; ...
 	NaN	2	2	1	3	3	NaN	1	3	NaN	NaN	1; ...
	NaN	2	2	1	3	3	NaN	1	3	NaN	NaN	1; ...
 	NaN	2	2	1	3	3	NaN	1	3	NaN	NaN	1; ...
 	NaN	2	2	NaN	NaN	NaN	NaN	1	3	NaN	NaN	1; ...
 	NaN	2	2	NaN	NaN	NaN	NaN	1	3	NaN	NaN	1; ...
 	NaN	2	2	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	1; ...
 	NaN	2	2	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	1];

%V	2	3	4	2	3	4	2	3	4	2	3	4
V=[ NaN	2	2	1	3	3	1	1	3	NaN	NaN	1; ...
    NaN	2	2	1	3	3	1	1	1	NaN	NaN	1; ...
    NaN	2	2	1	3	3	1	1	1	NaN	NaN	1; ...
    NaN	2	2	1	3	3	1	1	3	NaN	NaN	1; ...
    NaN	2	2	1	3	3	1	1	1	NaN	NaN	1; ...
    NaN	2	2	1	3	3	1	1	3	NaN	NaN	1; ...
    NaN	2	2	1	3	3	1	1	1	NaN	NaN	1; ...
    NaN	2	2	1	3	3	1	1	1	NaN	NaN	1; ...
    NaN	2	2	1	3	3	1	1	3	NaN	NaN	1; ...
    NaN	2	2	1	3	3	1	1	3	NaN	NaN	1; ...
    NaN	2	2	1	3	3	1	1	3	NaN	NaN	1; ...
    NaN	2	2	1	3	3	1	1	3	NaN	NaN	1; ...
    NaN	2	2	1	3	3	1	1	3	NaN	NaN	1; ...
    NaN	2	2	1	3	3	1	1	3	NaN	NaN	1; ...
    NaN	2	2	1	3	3	NaN	1	3	NaN	NaN	1; ...
    NaN	2	2	1	3	3	NaN	1	3	NaN	NaN	1; ...
    NaN	2	2	1	3	3	1	1	3	NaN	NaN	1; ...
    2	2	2	1	3	3	1	1	3	NaN	NaN	1; ...
    2	2	2	1	3	3	1	1	3	NaN	NaN	1; ...
    2	2	2	1	3	3	1	1	3	NaN	NaN	1; ...
    2	2	2	1	3	3	1	1	3	NaN	NaN	1; ...
    2	2	2	1	3	3	1	1	3	NaN	NaN	1; ...
    2	2	2	1	3	3	1	1	3	NaN	NaN	1; ...
    2	2	2	1	3	3	1	1	3	NaN	NaN	1; ...
    2	2	2	1	3	3	NaN	1	3	NaN	NaN	1; ...
    2	2	2	1	3	3	NaN	1	3	NaN	NaN	1; ...
    2	2	2	1	3	3	NaN	1	3	NaN	NaN	1; ...
    2	2	2	1	3	3	NaN	1	3	NaN	NaN	1];


U_reshape=cell(n_years,1);
for h=1:n_years
    U_reshape{h}=NaN(n_types, n_types);
    for i=1:n_types
        U_reshape{h}(i,2:n_types)=U(h,(i-1)*(n_types-1)+1:i*(n_types-1));
    end
end

V_reshape=cell(n_years,1);
for h=1:n_years
    V_reshape{h}=NaN(n_types, n_types);
    for i=1:n_types
        V_reshape{h}(2:n_types,i)=(V(h,(i-1)*(n_types-1)+1:i*(n_types-1))).';
    end
end

%%
fp='.../SpecB/output'; %specify path

D_31=zeros(n_years,2); %[min max] for each year
D_32=zeros(n_years,2); %[min max] for each year
D_43=zeros(n_years,2); %[min max] for each year

D_21=zeros(n_years,2); %[min max] for each year
D_41=zeros(n_years,2); %[min max] for each year
D_42=zeros(n_years,2); %[min max] for each year


for h=1:n_years
    
    IdSetM=cell(n_types,1);
    IdSetW=cell(n_types,1);
    D_interval=cell(n_types,n_types);
    
    for j=1:n_types
        
       filepath = fullfile( fp, ['IdSetM' num2str(all_years(h)) '_final_' num2str(j) '.mat'] );
       if exist( filepath, 'file' )
          A = getfield( load( filepath ), 'IdSetM_final' ); 
          IdSetM{j}=A;
       end 
       
       filepath = fullfile( fp, ['IdSetW' num2str(all_years(h)) '_final_' num2str(j) '.mat'] );
       if exist( filepath, 'file' )
          A = getfield( load( filepath ), 'IdSetW_final' ); 
          IdSetW{j}=A;
       end  
    end
    
    for i=1:n_types-1
    for k=i+1:n_types
        %Men
        M_ii_ik=unique([IdSetM{i}(:,i+1) IdSetM{i}(:,k+1)], 'rows', 'stable'); 
        M_kk_ki=unique([IdSetM{k}(:,k+1) IdSetM{k}(:,i+1)], 'rows', 'stable');
        [cx1, cx2] = ndgrid(1:size(M_ii_ik,1), 1:size(M_kk_ki,1));
        cx1=cx1(:);
        cx2=cx2(:);
        M_ii_kk_ik_ki=[M_ii_ik(cx1,1) M_kk_ki(cx2,1) M_ii_ik(cx1,2) M_kk_ki(cx2,2)]; %[U_ii U_kk U_ik U_ki]
        
        %Women
        W_ii_ki=unique([IdSetW{i}(:,i+1) IdSetW{i}(:,k+1)], 'rows', 'stable'); 
        W_kk_ik=unique([IdSetW{k}(:,k+1) IdSetW{k}(:,i+1)], 'rows', 'stable');
        [cx1, cx2] = ndgrid(1:size(W_ii_ki,1), 1:size(W_kk_ik,1));
        cx1=cx1(:);
        cx2=cx2(:);
        W_ii_kk_ik_ki=[W_ii_ki(cx1,1) W_kk_ik(cx2,1) W_kk_ik(cx2,2) W_ii_ki(cx1,2)]; %[V_ii V_kk V_ik V_ki]
        
        %Compose
        %A(h,1)+B(:,1)+A(h,2)+B(:,2)-(A(h,3)+B(:,3))-(A(h,4)+B(:,4))
        BB = W_ii_kk_ik_ki(:,1) + W_ii_kk_ik_ki(:,2) - W_ii_kk_ik_ki(:,3) - W_ii_kk_ik_ki(:,4);
        AA = M_ii_kk_ik_ki(:,1) + M_ii_kk_ik_ki(:,2) - M_ii_kk_ik_ki(:,3) - M_ii_kk_ik_ki(:,4); 
        D_interval{i,k}= [min(BB) max(BB)] + [min(AA) max(AA)];

        
        
        %Replace +Inf, -Inf
        if (U_reshape{h}(i,i)==1 || U_reshape{h}(i,i)==3)  || ...
           (U_reshape{h}(k,k)==1 || U_reshape{h}(k,k)==3)  || ...
           (V_reshape{h}(i,i)==1 || V_reshape{h}(i,i)==3)  || ...
           (V_reshape{h}(k,k)==1 || V_reshape{h}(k,k)==3)  || ...
           (U_reshape{h}(i,k)==2 || U_reshape{h}(i,k)==3)  || ...
           (U_reshape{h}(k,i)==2 || U_reshape{h}(k,i)==3)  || ...
           (V_reshape{h}(i,k)==2 || V_reshape{h}(i,k)==3)  || ...
           (V_reshape{h}(k,i)==2 || V_reshape{h}(k,i)==3)     
           D_interval{i,k}(2)=Inf;
        end
       if ( U_reshape{h}(i,i)==2 || U_reshape{h}(i,i)==3)  || ...
           (U_reshape{h}(k,k)==2 || U_reshape{h}(k,k)==3)  || ...
           (V_reshape{h}(i,i)==2 || V_reshape{h}(i,i)==3)  || ...
           (V_reshape{h}(k,k)==2 || V_reshape{h}(k,k)==3)  || ...
           (U_reshape{h}(i,k)==1 || U_reshape{h}(i,k)==3)  || ...
           (U_reshape{h}(k,i)==1 || U_reshape{h}(k,i)==3)  || ...
           (V_reshape{h}(i,k)==1 || V_reshape{h}(i,k)==3)  || ...
           (V_reshape{h}(k,i)==1 || V_reshape{h}(k,i)==3)  
           D_interval{i,k}(1)=-Inf;
       end
           
    end
    end
    
    D_31(h,1)=D_interval{1,3}(1);
    D_31(h,2)=D_interval{1,3}(2);
    D_32(h,1)=D_interval{2,3}(1);
    D_32(h,2)=D_interval{2,3}(2);
    D_43(h,1)=D_interval{3,4}(1);
    D_43(h,2)=D_interval{3,4}(2);
    
    D_21(h,1)=D_interval{1,2}(1);
    D_21(h,2)=D_interval{1,2}(2);
    D_41(h,1)=D_interval{1,4}(1);
    D_41(h,2)=D_interval{1,4}(2);
    D_42(h,1)=D_interval{2,4}(1);
    D_42(h,2)=D_interval{2,4}(2);
    
    
end


save('D_R2', 'D_31', 'D_32', 'D_43', 'D_21', 'D_41', 'D_42')

